%% Random Number Generator

%help rng

%Generated numbers are pseudo-random.  They are produced by an algorithm. If
%the algoritm is seeded the same way, the same numbers will be produced:

rng(0)
rand(3,1)

rng(0)
rand(3,1)

%This can be use to generate the same set of 'random' numbers each time.  
% See 'doc rng' for algorithm options. 


%% The Inversion Method

web('https://en.wikipedia.org/wiki/Inverse_transform_sampling')

%An example for the exponential distribution:
% PDF: lambda*exp(-lambda*x)
% CDF: 1-exp(-lambda*x)

N=1000;
lambda=1;

%This is the inverse of the CDF.  
Data=-log(1-rand(N,1))/lambda;

figure
hist(Data,20)
xlabel('X')
ylabel('Counts')
%Note the data looks to be correctly distributed. 

%The mean of the PDF should be lambda. Check:
fprintf('Mean should be %g.  Mean of Data is %g \n',lambda,mean(Data))
%The variance of the PDF should be lambda^-2. Check:
fprintf('Variance should be %g.  Variance of Data is %g \n',lambda^2,var(Data))


%We also briefly talked about rejection sampling:
web('https://en.wikipedia.org/wiki/Rejection_sampling')




%% The Poisson Distribution

%We get Poisson distribution when we are counting the number of events that
% occur in a fixed time period from a constant-rate process.  This is also
% called a 'Poisson' process. 

web('https://en.wikipedia.org/wiki/Poisson_distribution')

TSteps=1000;  %Time steps
P=.020;  %Probability of an event during time interval
NTrials=10000;

%Vary P from small to large to change expected value from ~ 1 to > 10

%Expected, or mean number of events over repeated trials
Expected=P*TSteps

%This is the events per time step 
%Look carefully at what this does.  '1' is an event. 
Events=double(rand(TSteps,NTrials)<P);

%total counts in each experiment
Counts=sum(Events)


%Show the events. 
figure
hist(Counts,(min(Counts):max(Counts)))
xlabel('Number of Events per Trial')
ylabel('Occurance')




